home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / answrbok / 4_10.lha / 4_10 / 4_10b.c < prev    next >
Text File  |  1993-08-08  |  1KB  |  66 lines
  1. * Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
  2. * The C++ Answer Book */
  3. * Tony Hansen */
  4. * All rights reserved. */
  5. ifndef VER2        /* DELETE */
  6. include "4_10b0.c"    /* DELETE */
  7. else            /* DELETE */
  8. include "4_10c.c"    /* DELETE */
  9. endif            /* DELETE */
  10. / invert A and place the result in I
  11. nt matinv(double **A, double **I, int nelem)
  12.  
  13.    // initialize identity matrix I
  14.    for (int i = 0; i < nelem; i++)
  15. {
  16. for (int j = 0; j < nelem; j++)
  17.     I[i][j] = 0.0;
  18.  
  19. I[i][i] = 1.0;
  20. }
  21.  
  22.    // loop diag from 0 to nelem-1
  23.    for (int diag = 0; diag < nelem; diag++)
  24. {
  25. if (!dopivot(A, I, diag, nelem))
  26.     return 0;
  27.  
  28. // divide row diag by A[diag][diag] in A and I
  29. //     ignore A[diag][0..diag-1] since they
  30. //     are already 0
  31. double div = A[diag][diag];
  32. if (div != 1.0)
  33.     {
  34.     A[diag][diag] = 1.0;
  35.     for (int j = diag+1; j < nelem; j++)
  36.     A[diag][j] /= div;
  37.  
  38.     for (j = 0; j < nelem; j++)
  39.     I[diag][j] /= div;
  40.     }
  41.  
  42. // loop i from 0 to diag-1 and
  43. // from diag+1 to nelem-1
  44. for (i = 0; i < nelem; i++)
  45.     {
  46.     if (i == diag)
  47.     continue;
  48.  
  49.     // subtract A[i][diag] * row diag
  50.     // from row i in A and I
  51.     double sub = A[i][diag];
  52.     if (sub != 0.0)
  53.     {
  54.     A[i][diag] = 0.0;
  55.     for (int j = diag+1; j < nelem; j++)
  56.         A[i][j] -= sub * A[diag][j];
  57.  
  58.     for (j = 0; j < nelem; j++)
  59.         I[i][j] -= sub * I[diag][j];
  60.     }
  61.     }
  62. }
  63.  
  64.    return 1;
  65.  
  66.